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We explore some consequences of the "alpha model," also called the "Lagrangian-averaged" model, for 
two-dimensional incompressible magnetohydrodynamic (MHD) turbulence. This model is an extension of the 
smoothing procedure in fluid dynamics which filters velocity fields locally while leaving their associated vor- 
ticities unsmoothed, and has proved useful for high Reynolds number turbulence computations. We consider 
several known effects (selective decay, dynamic alignment, inverse cascades, and the probability distribution 
functions of fluctuating turbulent quantities) in magnetofiuid turbulence and compare the results of numeri- 
cal solutions of the primitive MHD equations with their alpha-model counterparts' performance for the same 
flows, in regimes where available resolution is adequate to explore both. The hope is to justify the use of 
the alpha model in regimes that lie outside currently available resolution, as will be the case in particular in 
three-dimensional geometry or for magnetic Prandtl numbers differing significantly from unity. We focus our 
investigation, using direct numerical simulations with a standard and fully parallelized pseudo-spectral method 
and periodic boundary conditions in two space dimensions, on the role that such a modeling of the small scales 
using the Lagrangian-averaged framework plays in the large-scale dynamics of MHD turbulence. Several flows 
are examined, and for all of them one can conclude that the statistical properties of the large-scale spectra are 
recovered, whereas small-scale detailed phase information (such as e.g. the location of structures) is lost. 

PACS numbers: 47.27.Eq; 47.27.Gs; 47.11.+j 



I. INTRODUCTION 

One of the most persistent difficulties in the computation of 
the turbulent behavior of fluids and magnetofluids is the wide 
range of dynamically interacting length and time scales that 
have to be evaluated. At large Reynolds number, many orders 
of magnitude in length scales are implied, for example, in the 
dynamical behavior of the atmosphere, the oceans, or the so- 
lar wind, to take some familiar situations. For many purposes, 
it might be adequate to compute only the long-wavelength 
components of the spectra of the fields involved if some more 
economical representation or model of the small scale behav- 
ior could be given which would not do violence to the accu- 
racy with which the large scales are computed. Such topics 
as "large eddy simulation" and "eddy viscosity," designed to 
cope with this difficulty, have generated a vast literature, one 
which we make no attempt to survey here (see e.g. Refs. (Ji- 
ll)- 

A novel offering along these lines which has appeared in 
recent years is the so-called "alpha model" of Holm, Foias, 
Margolin, Marsden, Olson, Ratiu, Titi, Wynne and especially 
Chen, whose comparisons with turbulent channel and pipe 
flow called the most attention to the alpha model's possi- 
bilities {e.g. , Refs. 0|-i3l; many other references could 
also be cited). The model is also variously called the "La- 
grangian averaged model" or, in some of the earliest papers, 
the "Camassa-Holm" equations. This alpha model is subject 
to a variety of derivations, interpretations and connections, 
ranging from the mathematically sophisticated 0, Q3 HUl to 
the intuitive and simple 1 1311 . It is of interest to subject its 
predictions to tests against both experimental data (see Ref. 
|8[) and numerical solutions of the relevant continuum equa- 
tions to which alpha modeling has not yet been applied (see 
e.g. Ref. |21 for the Navier-Stokes equations in three dimen- 
sions). Its extension to the case of coupling to a magnetic 



field in the magnetohydrodynamic (MHD)limit and in the 
non-dissipative case can be found in 11411 . In that con- 
text, the main purpose of this article is to carry out some of 
the numerical tests that have not previously been done for the 
case of MHD. 

It is to be emphasized that there is no derivation of the al- 
pha model that is completely systematic and deductive. Every 
presentation of it has involved steps that call for justification 
by their consequences, and that is the spirit in which we are 
proceeding here. In Ref. 11311 . which seems the most eco- 
nomical derivation possible, the point of view is taken that 
we smooth the fields (e.g. , the velocity field v and, in MHD, 
the magnetic field B) but not their "sources" (e.g. , the vor- 
ticity field uj and, in MHD, the electric current density j). By 
"sources," we mean here the curls of v and B which, given a 
set of boundary conditions, determine them through the Biot- 
Savart law or solutions to Poisson's equation. The assumption 
is that the large-scale dynamics are relatively insensitive to 
small changes in the positions of the sources, but are sensi- 
tive to the strengths and approximate positions of them. Said 
another way, the large-scale fields alone are assumed to be 
responsible for those motions of their sources which signifi- 
cantly affect those large-scale fields. This assumption is by no 
means self-evident, but does have the advantage of reducing 
the derivation to a single algebraic step, in contrast to some 
more involved derivations which have been given, and which 
seem logically no more compelling. Our focus here is on nei- 
ther presenting unarguable alpha-model derivations or com- 
paring the possible variants of it, but rather on exploring the 
consequences of the primitive version of it given here [Eqs. 
il 11 and d!2i in what follows]. 

Modeling MHD flows with the Lagrangian-averaged 
methodology has been barely explored with no emphasis on 
the turbulence regime. This article thus focuses on several 
such predictions, numerically obtained, for the case of two- 
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dimensional magnetohydrodynamics, or hereafter 2D MHD 
(see e.g. for a brief review, Ref. I15l0 . Several effects have 
been studied phenomenologically, theoretically and numeri- 
cally in the past, and may be identified in the literature by 
the names selective decay, dynamic alignment, direct and in- 
verse cascades, and the characterization of probability dis- 
tribution functions (pdfs) for the fluctuating field variables. 
There are several Reynolds-like numbers which can be at- 
tributed to MHD turbulent flows, since there are two possible 
velocities which may appear in the numerators (the flow speed 
and the Alfven speed) and two diffusivities that may appear 
in the denominators (kinematic viscosity and magnetic diffu- 
sivity). Some length scale characteristic of the initial fields 
is usually present in the numerators. All these Reynolds-like 
numbers can be made to appear in the places of reciprocals 
of the transport coefficients in front of the dissipative terms in 
various dimensionless representations of the MHD equations. 
In general, the larger the values of these Reynolds-like num- 
bers (or equivalently, the smaller the transport coefficients), 
the greater the required numerical resolution to follow their 
solutions. Values of a Reynolds number like 10 4 usually strain 
available computer resources even in two space dimensions 
(2D), and while the attainable total number of degrees of free- 
dom with computers has been steadily increasing over several 
decades, there are situations in which one might be curious 
about results in cases of far higher values of direct interest 
for geophysical flows and yet not attainable in the foresee- 
able future. The alpha model, if it can be verified to give cor- 
rect predictions in the range of accurate, un-modelized solu- 
tions, will acquire a certain credibility in providing the behav- 
ior (at least of the long-wavelength Fourier components) in 
situations with Reynolds-like numbers so high as to put them 
presently far out of reach of direct numerical solutions (DNS), 
particularly so in three space dimensions (3D). Another set of 
regimes where modeling is needed is when widely disparate 
time and length scales occur, such as for either a small or large 
magnetic Prandtl number; this is the case for the former in liq- 
uid metals as encountered in laboratory dynamo experiments 
and in breeder reactors, in the core of the earth and planets, 
and in the convective zones of the sun and stars, or for the 
latter in the interstellar medium. 

In Section II, we write down the alpha-model equations that 
we assume for incompressible, one-fluid MHD with a min- 
imum of theoretical justification, following Ref. |13]; they 
include the effects of true viscous and dissipation. We pro- 
vide expressions for ideal invariants that are conserved by the 
alpha model when the viscous and Ohmic dissipation coeffi- 
cients are dropped, and decay laws for them when the dissipa- 
tion coefficients are present and finite. Much of the approach 
to and vocabulary of the way turbulence problems have his- 
torically been formalized for MHD (a subject where computa- 
tional data vastly exceeds experimental or observational data) 
can by now be taken for granted. In Section EH, we describe 
results for the problem of selective decay in 2D. In Section 
IV, we turn our attention to that of dynamic alignment of the 
velocity and magnetic fields in turbulent decays. In Section V, 
we focus on inverse cascade computations. In all these cases, 
there are comparisons to be made between the alpha-modeled 
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FIG. 1 : (a) Magnetic energy (upper curves) and kinetic energy (lower 
curves) as a function of time until t = 5, and (b) cross helicity as a 
function of time until t — 20, for the selective decay runs 1-4 (see 
Table . The temporal evolutions of DNS and alpha runs appear 
similar. 



results and direct solutions of the primitive MHD equations. 
In Section VI, we address ourselves to the problem of quanti- 
tatively assigning errors to the alpha model, as compared with 
well-resolved full MHD as well as with unresolved MHD, of 
resolution comparable to that used for the alpha model. In 
Section VII, we sample a few effects at Reynolds-like num- 
bers that are too high for any immediately foreseeable DNS 
code to approach. Finally, in Section VIII, we briefly summa- 
rize the results and suggest future problems in which the alpha 
model may have some utility. 

To anticipate the conclusions of the paper, we note some 
features of the DNS solutions that the alpha model apparently 
finds out of reach. For instance, the location, in the plane, of 
specific features of evolving turbulent fields such as contour 
plots of vorticity or vector potential are virtually never accu- 
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rately reproduced after short times (contours of constant mag- 
netic vector potential in 2D are magnetic field lines). Like- 
wise, the spectral details at the small scales are not accurate 
and are not expected to be, since it is modifications of the dy- 
namics at small scales that make the alpha model possible in 
the first place. But as far as the long wavelength component 
behavior for the turbulent kinetic and magnetic spectra is con- 
cerned, the alpha model seems to recover the main features of 
MHD turbulent flows in two space dimensions. 

One technical feature of the computations peculiar to two- 
dimensional (2D) MHD should be commented upon. Because 
of the inherent tendency of 2D magnetic excitations to mi- 
grate to longer wavelengths, treatment of initial value prob- 
lems requires beginning with excitations located in interme- 
diate length scales, rather than at the longest wavelengths. In 
wavenumber space, this means that any filtering that is done 
must be done above wavenumbers corresponding to shorter 
wavelengths than those in the initial conditions. For this rea- 
son, very low resolution alpha-model calculations are not pos- 
sible, unless the initial conditions themselves are to be left 
outside the basic box in Fourier space. This contrasts with 
the situation in three-dimensional hydrodynamics, where the 
emphasis is typically on cascades to shorter length scales, and 
where for the above reasons, it is possible to attempt large 
eddy simulations (LES) with very low maximum wavenum- 
bers, and the initial excitations may all reside at the largest 
scales. While the ratio of our maximum retained wavenum- 
ber to the wavenumber where the filtering begins is about 
8, it should be noted that larger ratios are feasible for three- 
dimensional Navier-Stokes LES and alpha model computa- 
tions. 



II. THE ALPHA MODEL FOR MHD 

A. The equations 

We write the three dimensional version of the equations first 
for the primitive incompressible MHD equations, and then for 
the alpha model. We will then specialize them to two dimen- 
sions for the purposes of this paper. The basic variables are the 
velocity field v and the magnetic field B, functions of space 
and time coordinates (x, t). In dimensionless (Alfvenic) units, 
the equations are: 



<9v 

~dt 
dB 

~dt 



+ v- Vv = -VP+jxB-i/Vxu (1) 



v • VB = B • Vv - 77V x j 



(2) 



together with V • v = and V • B = 0. 

The velocity field may be considered to be expressed in 
units of an r.m.s. value of the initial fluctuating velocity field, 
which we typically take to be 1 . The magnetic field is made di- 
mensionless by solving for the magnetic field value that would 
lead to an Alfven speed equal to the r.m.s. velocity field and 
dividing the magnetic field in laboratory units by that. The 
mechanical pressure is V, which has first been divided by the 



mass density and then expressed in the units of the dimension- 
less velocity. The mass density is assumed to be constant and 
uniform. The viscosity v and the magnetic diffusivity 7/ can be 
considered to be reciprocals of the mechanical and magnetic 
Reynolds numbers, respectively, in these units. Anticipating 
that the computations will be carried out inside a periodic box 
of edge 2tt, the unit of length will in general be taken to be 
equal to unity, or about 1/6 of a box dimension. 

In the dimensionless units, the curl of the velocity is u>, 
the vorticity field, and the curl of the magnetic field is j, the 
electric current density. The magnetic field B can be written 
as the curl of a vector potential A, which, removing a curl 
from Eq. obeys 



dA 
~dt 



= v x B — 77J — V$ 



(3) 



where the scalar potential is $. $ can be determined by taking 
the divergence of Eq. (|3}, imposing the Coulomb gauge on A 
(i.e. writing V • A — 0), and solving the resulting Poisson 
equation for $, involving v and B in the source term. (In a 
similar way, the pressure V can be found by taking the diver- 
gence of Eq. 0, using the vanishing of the divergence of the 
time derivative of the velocity field v, and solving the result- 
ing Poisson equation for the pressure. These Poisson solutions 
are easy to solve in Fourier space.) 

To obtain the alpha-model version of Eqs. ( 1 1 131 . we may 
divide v and B into smoothed values plus fluctuations about 
those values. Thus 



and 



v = u, + 5v 



B = B s + <5B 



(4) 



(5) 



Here, the smoothed values of the fields, u s and B s , are 
defined by 



, ,0 , expl — |x — x'l/al , , . 
MX Wlx-x V(X '^ (6) 



B s = d 



c , exp[-|x - x'|/a] 
47ra 2 |x — x'l 



B(x',i) 



(7) 



with a at this point an arbitrary length; a -1 is typically to be 
chosen as larger than the wavenumbers whose behavior it is 
desired to reproduce accurately. 

In general, we use the same values of a for smoothing v 
and B, though we note the possibility of assigning a priori 
one value for v and a different value for B (respectively, ctk 
and a m ). Another possibility is to choose a m — for B, 
which in that case will leave us with an unsmoothed magnetic 
field. We will show such an example later. 

We now take the curl of Eq. Q to obtain the equation of 
motion in the vorticity representation, 



~dt 



v-Vu = u- Vv + Vx(jxB) + isV 2 u) (8) 
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FIG. 2: (a) Magnetic energy and (b) kinetic energy spectra, for selective decay runs 1-4 (see Table[I}, at t = 5; and (c) magnetic energy 
and (d) kinetic energy spectra at t = 100. The vertical line gives k a ~ 1/a. The crosses (+) indicate the place on the spectrum where an 
under-resolved DNS run (256 2 grid points) departs significantly from the resolved computed DNS spectra. This convention will also be used 
in subsequent figures: crosses always indicate part of the under-resolved DNS spectrum for the same initial conditions and times. 
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and then substitute into Eqs. Q and l|8} the fields expressed 
in Eqs. @ and (|5j- Note that we do not smooth the vorticity 
U3 which can be regarded as the source, in a Poisson or Biot- 
Savart sense, of v. Nor do we smooth j, which bears the same 
mathematical relation to B as w does to v. The result is: 



FIG. 3: Square vector potential energy spectrum shown at t = 300, 
for the same selective decay runs as in Fig. [2] By this late time, 
essentially all of A is concentrated in k = 1. 
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did n 

— + (u s + <5v) • Vw = u) ■ V(u s + 8v) + V x [j x (B s + SB)] + isV a uj 



(9) 



and 



9 t (B s + <5B) + (u s + 5v) ■ V(B S + SB) = (B s + SB) ■ V(u s + 5v) - 77V x j 

I 



(10) 



upon which no approximations have as yet been made. That 
is, they are equivalent to Eqs. Q-(|2ji. 

Taking a modeling or heuristic point of view fl3ll . the 
essence of the alpha model is to neglect the fluctuations Sv 
and SB in relation to the smoothed fields u s and B s in Eqs. 
(|9) and dl0> . while leaving the source terms u> and j alone. 
This is one way of looking at the alpha approximation. Its re- 
lation to other, more complicated derivations will not be dis- 
cussed here, since our intent is to test the alpha model rather 
than to justify it from anything like first principles. Further 
discussion of the above approximation, which is the only one 
in our formulation, can be expected elsewhere. 

The alpha model equations are then (removing a curl from 
Eq. ©): 



Noting that only one component of A, the z-component, is 
relevant to two dimensions, the result is: 



■^ + u s -Vw = B s - Vj + i/V 2 w, (15) 
at 



dA Sk 
dt 



u s • VA H . = 



-m » 



(16) 



where there are stream functions ip and vector potentials 
A z that bear Poisson relations to their sources, both for the 
smoothed and unsmoothed versions: 



<9 f v + u s • Vv+VjVu J a = -VP ljxB s -vVxu (11) 
and 

d t B s + u s • VB S = B s • Vu s -ijVxj. (12) 

Note that the smoothed quantities bear the subscript letter s, 
and the unsmoothed ones do not. We shall follow this con- 
vention throughout. Eq. (II 21 could be viewed alternatively 
as a hyper-resistivity approximation on B s . The connection 
between the smoothed and unsmoothed fields may be stated 
in differential form as 

v = (l-a 2 V 2 )u s (13) 

and 

B = (1 -a 2 V 2 ) B s . (14) 

We may associate smoothed values of u> and j with the un- 
smoothed ones according to the same recipe; even though 
they do not enter directly into the dynamical equations, they 
are at some points convenient to think and talk about. Thus 
V x u s = u> s , similarly V x B s = j s , and V x A s = B s . 
A smoothed vector potential A s may be regarded as having 
a curl B s , while obeying a Poisson relation to j s , namely 
V 2 A S = — j s . We stress that ui s and j s do not enter the alpha 
model equations we use. 

Specialization to two dimensions is achieved by taking the 
curl of Eq. Q or specializing Eq. and Eq. (0 to the 
two dimensional geometry in which there are only two (x, y) 
non-zero components of v and B, and only one component 
(z) of us or j, and carrying out the smoothing approximations 
so described. All fields are independent of the z coordinate. 



V 2 * = -u , V 2 * s = -lu s , (17) 



W 2 A Z = -j , V 2 A Sz = -j s . (18) 

To re-iterate, the principal intent of this paper is to compare 
typical solutions of Eqs. d!5l > and with solutions, for the 
same initial and boundary conditions, of the well-known 2D 
MHD equations, unsmoothed. 

B. The invariants 

There are three ideal invariants for both sets of equations: 
the total energy E, the total cross helicity He, and the total 
mean-square vector potential A. The alpha model expressions 
for these are, respectively, 



E=- / d 2 x(u s -v + B s -B), (19) 



H C = 2 j rf 2 xv -B s , (20) 

and 



A=-jd 2 xAl. (21) 

Note that the energy invariant E involves both the smoothed 
and unsmoothed velocity and magnetic field, whereas only the 
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TABLE I: Main characteristics of the runs, and a^ 1 are the re- 
ciprocal of the alpha lengths for the velocity and the magnetic field; 
TV is the grid resolution before dealiasing, R\ is the Taylor Reynolds 
number (see Eq. 1251 at peak dissipation, and the last column gives 
the figures relating to the different runs, namely: runs 1-4 for selec- 
tive decay, runs 5-7 for dynamic alignment, runs 8-12 for the inverse 
cascade of magnetic potential and runs 13-17 for large-scale turbu- 
lence. In the figures, solid lines are for fully resolved DNS (runs 1, 
5, 8 & 13), dashed lines for runs 2, 6, 9 & 14, dashed-triple dots for 
runs 3, 11 & 15, dotted lines for runs 4, 7, 12 & 16, and a dash-dot 
line for run 10. 



smoothed magnetic variables appear in the expressions of He 
and A, due to the linearity of the induction equation, once 
the velocity field is given; however, the decay rates of He 
and A involve both the smoothed and unsmoothed magnetic 
variables whereas the decay rate of energy only involves the 
unsmoothed current density (see below). 

The decay laws for these quantities are, in periodic bound- 
ary conditions, 



dE 
~dt 



dH c 



dt 



Xujj - -v 



j2 -1 



d 2 x ujj s , 



(22) 



(23) 



and 

dA 
It 



r\ \ d x A S J 



->1 



d 2 x VA. • VA Z . (24) 



These three invariants will be at the core of the numeri- 
cal tests to be reported in the next three Sections. Finally, it 
should be noted that Eqs. l[T9}-(|24} differ from their full MHD 
equivalents in detail, but approach them as a — •» 0. 

For convenience and later reference, the main characteris- 
tics of all the runs described in this paper are given in Table |U 
together with the number of the figures related to the different 
category of runs. N is the grid resolution before dealiasing, 
using the standard 2 /3 rule in all the runs described in this pa- 
per (hence, for a run on a grid of N x N points, the maximum 
wavenumber attainable is equal to N/3). Finally, R\ is the 
Taylor Reynolds number defined as 



R\ — 



it is based on the r.m.s. velocity and on the Taylor scale 



A = 27rv / < v 2 > I < LU 2 > 



(25) 



(26) 



computed at the peak of the dissipation. The viscosity is 
equal to 5 x 10~ 4 for the selective decay runs 1-4, with ini- 
tial conditions with non-vanishing Fourier coefficients (see 
Eq. ( I27i ~) for wavenumbers between ki = 10 and k^ = 30. 
Runs 5-7 are for dynamic alignment, with v = and 
ki = 5, &2 = 10. Runs 8-12 are for the inverse cascade of 
magnetic potential, with v = 10~ 3 , the forcing occurring in 
the interval k\ = 18, = 22. Runs 13-16 deal with large- 
scale turbulence with i/ = 5x 10~ 4 and the initial conditions 
confined between k\ = 1 and &2 = 3, whereas for run 17, 
i/ = 2x 10 -5 . In the figures, solid lines are for fully resolved 
DNS (runs 1, 5, 8 & 13), dashed lines for runs 2, 6, 9 & 14, 
dashed-triple dots for runs with ak 7^ 0, a m = (runs 3,11 
& 15); finally, dotted lines are for runs 4, 7, 12 & 16, and a 
dash-dot line for run 10. Note that all runs have unit magnetic 
Prandtl number; the initial conditions are such that kinetic and 
magnetic energies are equal and of order unity with random 
phases. All runs are decaying {i.e. no forcing), except for the 
inverse cascade runs, which have zero initial conditions and 
forcing in the induction equation only. 



III. SELECTIVE DECAY 

By "selective decay," we mean turbulent processes (see e.g. 
1 16]-| 18]) in which one or more ideal invariants are dissipated 
rapidly relative to another, due to the transfer of the dissipated 
quantities to short wavelengths where the dissipation coeffi- 
cients become effective. In 2D MHD, with negligible cross- 
helicity, the selectively dissipated quantity is energy, while 
the nearly-conserved quantity is mean square vector potential. 
The limit defines a variational problem which seeks the state 
in which the dissipated quantities are as close to zero as they 
can be for the surviving value of the nearly-conserved quan- 
tity. There are no constraints on the cascade of kinetic energy 
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FIG. 4: Contour plots of vector potential (left) and current density (right) at t — 400 for the DNS run (top) and alpha run 3 (middle) and 4 
(bottom), with positive and negative values respectively in solid and dashed lines. Large scales are almost in the form of bars parallel to the 
axes. 



to short wavelengths, so the asymptotic state is expected to be 
one that is largely magnetic and has the surviving magnetic 
excitations peaked at the longest wavelengths. In particular, 
the vector potential spectrum should have a sharp maximum 
at the lowest wavenumber of the computation, here k m i n = 1. 
This effect has been demonstrated repeatedly in the past (| 16]- 
E). 

In this Section, we compare the full MHD (a = 0) results 
for a selective decay run of a familiar type with the conse- 



quences of the alpha model for the same initial conditions but 
finite a. We specify the initial conditions in Fourier space, 
with v and B represented as the Fourier series, 



v(x, t)=J2 v(k, t)e tk - x , B(x, t) = B(k, t)t 



ik-x 



(27) 

with similar decompositions for the other vector fields. The 
non-vanishing initial Fourier coefficients are confined to a ring 
in k-space between fci = 10 and k-i — 30. The amplitudes in 
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FIG. 5: Pdfs at t = 15 of (a) the current density, and (b) ln(^oj 2 + 
r]j 2 ), for all selective decay runs (see Table[I}. Note the larger values 
of current density for the alpha runs compared to the DNS (solid 
line). 



this ring are chosen equal, for both v and B, and the phases 
are chosen from a random-number generator. The overall nor- 
malization is such that the initial kinetic and magnetic ener- 
gies are both 0.5, referred to a unit (two-dimensional) volume. 
The maximum value of k is, after the de-aliasing, k max = 341 
and the time step is At = 5 x 10~ 4 ; the dimensionless vis- 
cosity and resistivity are both equal to 5 x 10~ 4 ; the initial 
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FIG. 6: (a) Temporal evolution of magnetic energy (top curves) and 
kinetic energy (bottom curves), and (b) of normalized cross helicity, 
for all dynamic alignment runs (see Table[I). Note that only the first 
five units of time for the run are shown in (a); the solid line is for the 
DNS run. 



Reynolds numbers, kinetic and magnetic, are then formally 
equal respectively to 2000, based on a unit length scale in the 
basic box of edge 27r, whereas the Taylor Reynolds number 
at peak of dissipation is equal to 215 for the DNS run, and 
slightly larger for the alpha runs (see Table [j}. The time will 
be measured in units that are defined by the ratio of unit length 
to the initial r.m.s. velocity; based on the energy containing 
scale, one time unit can be several initial eddy turnover times, 
a number which may increase or decrease as the kinetic en- 
ergy is dissipated. 

Figure ^ displays the computed magnetic energy (upper 
curves) and kinetic energy (lower curves) versus time, show- 
ing the ultimate dominance of the magnetic energy over the 
kinetic, by an order of magnitude at this time. 

The solid lines are the results of the full MHD computa- 
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FIG. 7: (a) Magnetic and (b) kinetic energy spectra, for the three dynamic alignment runs (see Table[I}, at t = 2; (c) Magnetic and (d) kinetic 
energy spectra at t = 30. As elsewhere, crosses (+) indicate points on the computed under-resolved DNS spectrum (128 2 grid points) for the 
same initial conditions, and indicate the values of k where the first significant departures occur from the well-resolved DNS. 



tion (i.e. , a = 0), with 1024 2 grid points. The dashed lines 
are the results of the alpha model with 1024 2 grid points and 
with a = 1/40. The dotted line (barely distinguishable from 
the dashed one) shows the results for an alpha-model run with 
256 2 grid points and the same value of a. Note that a full 
MHD run with the lower number of grid points, i.e. an under- 
resolved computation of MHD turbulence, would display dis- 
agreement with the other three runs. We shall discuss this 
question in Section VI. 

Fig. \\]p shows the cross helicity (which remains very small, 
relative to the energy, because the random phases for the two 
fields imply negligible correlation between them) as a func- 
tion of time. The alpha model computations disagree with the 
exact MHD run, but since all the quantities are so small, this 
disagreement is not deemed to be significant, but rather oc- 
curs as fluctuations around values close to zero. We should 
remark that in both cases, the comparison between the full 
MHD quantities and their alpha-model analogues has been 
done after a rescaling of initial data which makes the ener- 
gies agree exactly at t = (with the full, unsmoothed MHD 
values); the original MHD energy is not quite the same as the 
energy integrals defined in Section II involving the smoothed 
fields, because of the smoothing, though the difference is only 
a few percent. 

Figures |2^,b show the omni-directional energy spectra for 



the magnetic and kinetic energies (as defined in Section II 
when a is non-zero), respectively, at t = 100. The same con- 
ventions adopted in Fig. [1] (i.e., solid lines mean full MHD, 
dashed lines mean alpha model with a = 1 /40 and 1024 2 grid 
points, and dotted lines mean runs done with the same value 
of a but with 256 2 grid points) will be followed throughout. 
The alpha-model spectra that are plotted result from taking 
the spectral density of the invariants in Equations dl9t- J21l i. 
The vertical line indicates the wavenumber k a corresponding 
to the length a. Note that the under-resolved spectra begin to 
differ at k ~ k a /2, but the a-model and well-resolved DNS 
agree up to k ~ k a . 

Figure |3 displays a spectrum in log-lin scales, of the vec- 
tor potential at very late times, when the selective decay is 
nearly complete and the magnetic excitations are concentrated 
in the longest wavelength allowed by the boundary conditions 
(kmin — 1)- The alpha model has reproduced this feature, 
with only a small disagreement in the values at k = 1 (see 
also Section VI for a more complete discussion of errors). 

The suppression of the small scales is quite apparent for 
both quantities, but the large scales, like the global energies 
exhibited in Fig.^, do not appear to be significantly affected. 
The lower resolution alpha model run reproduces the same 
result. This is what can be realistically hoped for from the 
alpha model, although we note that the disagreement between 
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FIG. 8: Right column: Contour plots of the (unsmoothed) stream function at t — 60 for dynamic alignment, with the DNS run at the top, run 6 
in the middle and run 7 at the bottom (see Table[lJ. Left column: contours of the (smoothed) vector potential at the same time. In all three runs, 
contours of magnetic potential and stream function are similar by that time, more so for the DNS run, and concentrated in the large scales. 



the true MHD run and the alpha runs starts at a scale roughly 
twice as large as oT 1 . 

Figures0]display contour plots of curves of constant vector 
potential A Sr (Fig. |4] left) and constant current density j (Fig. 
01 right) at time t = 100. The top panel is for full MHD, the 
middle one for a = 1/40 and 512 2 grid points, and the bot- 
tom one for a — 1/40 and 256 2 grid points. The flow may 
evolve toward a state reminiscent of those found in fl9ll in 
the case of 2D Navier-Stokes turbulence, with structures par- 



allel to either axis. While there are marked similarities in the 
kinds of structures present in the DNS and in the alpha runs, 
there are clearly no one-to-one correspondences as to specific 
features, either as to location, orientation, or intensity. From 
these and many similar figures we have looked at, we have 
concluded that while the alpha model does an excellent job of 
reproducing long-wavelength spectra, the pointwise details of 
the solution are not well tracked by it, at least in the absence 
of constraining material boundaries. 
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FIG. 9: Pdfs at t = 50 of (a) the current density, and (b) ln(^oj 2 + 
r)j 2 ), for all dynamic alignment runs (see Table 0. Note again the 
high values of |j| for the alpha runs. 



In Figs. 13 we display normalized probability distribution 
functions (pdfs) of the current density j in Fig. [5^, and in Fig. 
|5p of the spatial density of the dissipation rate of energy given 
in Eq. M2Y . note that the local (spatial) dissipation of kinetic 
energy differs in its expression, involving the symmetrized ve- 
locity gradient instead of the local squared vorticity density. 
The conventions with the lines are the same as those in the 
preceding three figures. It is apparent that the pdfs of the al- 



FIG. 10: (a) Temporal evolution of the magnetic energy (top curves) 
and kinetic energy (bottom curves) until t = 50, and (b) of the 
squared vector potential until t = 400, for all inverse cascade runs 
(see Table runs 8-12). Whereas energies are in good agreement 
with the DNS (shown as usual with a solid line), the growth of 
squared magnetic potential is slower for all alpha runs (see also Fig. 

El. 



pha model do a good job for the lower values of |j | but do 
not reproduce the tails accurately, in particular at lower res- 
olution, i.e. intermittency is not fully reproduced, and is not 
expected to be (for a study of intermittency in the context of 
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LES, see e.g. Ref. |20J]). The same is true of the dissipation 
density, although discrepancies appear smaller. 



IV. DYNAMIC ALIGNMENT 

A perfectly "aligned" solution to the ideal version of Eqs. 
Q and @ results whenever v = +B or v = B everywhere. 
Previous computations j2lll - ll23ll . inspired by observations in 
the quiet solar wind [24], have shown that MHD turbulence in 
which a significant degree of initial alignment, or correlation 
between the v and B fields, exists will evolve toward a state 
of greater and greater alignment as time goes on. The physi- 
cal origins of this process are not completely clear, except that 
we may note that an aligned state involves no spectral transfer 
to higher wave numbers where viscous and Ohmic dissipation 
are effective, so that those patches where alignment exists ini- 
tially may have a tendency simply to outlive the more active, 
unaligned patches where spectral transfer makes dissipation 
more likely. Similar alignment, this time between velocity 
and vorticity, can be observed for three-dimensional Navier- 
Stokes flows (see e.g. Ref. 12511 for an experimental study). 

A useful index of the global degree of alignment of v and B 
may be taken as the "fractional alignment": it is defined taking 
2Hc and dividing it by the square root of < u s • v >< B s • 
B >, where angle brackets mean spatial averages. When this 
ratio is unity, the fields may be regarded as perfectly aligned. 
In Fig. [6] we display the results of a run which starts with 
the Fourier amplitudes chosen to be equal in a ring with 5 < 
k < 10, with unit r.m.s. values of v and B and with phases 
chosen so that the fractional alignment is initially 0.3. Finally, 
At = 10~ 3 and v = rj = 10~ 3 ; the Reynolds numbers are 
equal to 1000, based upon unit length, unit r.m.s. velocity at 
t = and the transport coefficient, and the Taylor Reynolds 
number at peak dissipation is equal to 280 for the DNS run, 
and again slightly higher for the alpha runs (see Table[Q. 

In Fig. [6^, we show that the energies, magnetic (top curves) 
and kinetic (lower curves), as functions of time, have compa- 
rable evolutions. The fields, however, become progressively 
more aligned during their decay as can be seen in Fig. |6p, 
where the alignment index gradually increases from 0.3 to 
about 0.85. As before, the a = 0, or full MHD (with 512 2 grid 
points), results are exhibited as solid lines, the dashed lines are 
the results for a = 1/20 and 512 2 grid points, while the dot- 
ted line is for the same alpha but with only 128 2 grid points. 
The rather large discrepancies for long times for the corre- 
lation coefficient may come from the fact that it involves a 
normalization; the cross-helicity He themselves do not differ 
significantly (not shown); it also may imply that small scales, 
which are modified by the alpha modeling process, play a role 
in the growth of large-scale correlations between the velocity 
and the magnetic field. 

The spectra in Figs. 0show that in the aligning situation, 
the alpha model continues to do a good job of reproducing 
the long wavelength components of the MHD spectrum. Fig- 
ure|8]shows contour plots of the unsmoothed stream function 
(right) at t = 60, for full MHD (top), for the alpha model with 
a = 1/20 and 512 2 grid points (middle), and for the alpha 
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FIG. 11: (a) Squared vector potential, (b) magnetic, and (c) kinetic 
energy spectra, for all inverse cascade runs (see TableQ runs 8-12) 
at t — 160, with as usual a solid line for DNS runs; the two ver- 
tical arrows correspond to the two values of alpha. The line with a 
— 7/3 slope follows the phenomenological prediction for the mag- 
netic potential spectrum derived in 1 34]. Asa gets closer to the forc- 
ing wavenumbers, the inverse cascade is slowed down more signifi- 
cantly. 



model with a = 1/20 and 128 2 grid points (bottom); Figs. [8] 
(left) show the corresponding contours for the smoothed vec- 
tor potential. The pointwise alignment is more visible in the 
top contours (i.e. for the full MHD run) than in the finite a 
runs. 

Figure [9] shows pdfs of (a) the current density j, and (b) 
ln(r/j 2 + vuj 2 ) i.e. the spatial density of the energy dissipation 
rate; the units scales are lin-log; note the exponential depen- 
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FIG. 12: Flux of squared vector potential in Fourier space TlA(k) 
denned in Eq. J30l at t = 400: same runs and symbols as in Fig. 
II II Note the lesser amount of flux in the alpha runs compared to the 
DNS (solid line). 



dency at low values, as already found in 1 26], and also visible 
in Fig. [5] As in Section III, the only significant disagree- 
ment occurs in the tails, and the statistical properties of the 
turbulence are once again rather well reproduced by the alpha 
model approximation, even if the pointwise features shown in 
Fig. |8]are not particularly well reproduced. 



V. INVERSE CASCADE OF VECTOR POTENTIAL 

One might guess that one of the more demanding tests of 
alpha-modeled 2D MHD would be its performance in an in- 
verse cascade situation. The essence of the alpha model is 
that it suppresses the small scales to some degree. In direct 
cascade situations, the transfer is largely from the large scales 
to the small, and intuitively, it has often been reasoned that as 
long as the small scales can be made simply to disappear at 
that point, there should be no harm done in the large scales. 
Accurate or not, this reasoning is behind the use of Large 
Eddy Simulations (see e.g. Ref. |4] for a recent review), eddy 
viscosities and eddy resistivities that are sometimes employed 
to decrease the amount of resolution required for turbulence 
computations. In inverse cascades, it is the small scales which 
feed some cascadable quantity to the large scales, and any de- 
parture from proper dynamics at the small scales might be ex- 
pected to have non-trivial implications for the evolution of the 
large scales. 

The way inverse cascades have been studied numerically 
in the past, starting apparently with Lilly 12711 . is to write, 
on the right hand sides of the dynamical equations, forcing 
terms which are external to the fluid or magnetofluid equa- 
tions, and which are there to inject excitations in some field 
or another. These terms are typically band-limited in Fourier 
space, so that only a narrow range of fc-values, well above 



the energy-containing scales, are considered as externally ex- 
cited or stirred. The excitations can be injected either into the 
mechanical or the magnetic part of the dynamics. In stud- 
ies of the three-dimensional dynamo problem, the injection is 
typically into the velocity field, and involves the conversion 
of mechanical helicity into magnetic helicity, which is then 
transferred back into the long-wavelength part of the spectrum 
(see, e.g. , Refs. jH-jH). 

In two dimensions, the helicities are identically zero, and 
the anti-dynamo theorem prohibits the generation of persis- 
tent magnetic excitations by mechanical means, so the band- 
limited injections are magnetic and typically are considered 
to be the addition of mean square vector potential or magnetic 
flux, added randomly at the small scales (see, e.g. , Refs. IB2I1 - 
ll37lo . Here, we write a random forcing function on the right 
hand side of Eq. (1161 . which may be described as follows. We 
adopt a random forcing / only in the induction equation for 
the vector potential, of the form 



/(x,f)=^/(M)e 



ik-x 



(28) 



The sum runs from ki = 18 to k2 = 22. The amplitudes of 
all the coefficients /(k, t) in this ring are chosen equal, but the 
phases at each k are changed randomly with a correlation time 
larger than the time step, but smaller than the eddy turnover 
time. The phases are uniformly distributed between — 27r and 
27T. For the runs we will discuss in this Section, the correlation 
time was r = 2 x 1Q _1 , One would expect that the relation 
of the forcing band of wave numbers to the reciprocal of a 
would be a sensitive one in the outcome of an inverse vector 
potential cascade computation. This proves to be the case, 
and only in the situation where the forcing band lies at lower 
wave numbers than the reciprocal of a are recognizable results 
achieved. Even there, as will be seen in what follows, the 
agreement is less satisfactory than it has been for the selective 
decay and dynamic alignment tests. 

Figures UOh.b show the time histories of a magnetically 
forced run that started from an otherwise empty spectrum, 
with a time step At = 2 x 10~ 2 , and r; = v = 10~ 3 ; the 
Taylor Reynolds numbers at peak dissipation is for all runs 
~ 30. The upper curves in Fig. 110b are magnetic energies 
and the lower set are kinetic energies. The curves in Fig. El? 
are mean square vector potentials as functions of time. The 
solid lines are for a full MHD run (a = 0) with 256 2 grid 
points. The dashed lines are for a — 1/20 and 256 2 grid 
points. The dashed-dotted lines are for a = 1/30 and 256 2 
grid points. The dashed-triple-dotted lines are for the me- 
chanical afc = 1/20, but with the alpha parameter appearing 
in the induction equation a m set equal to zero (i.e. the mag- 
netic variables are unsmoothed). The dotted lines are for 128 2 
grid points and both alphas = 1 /30. The forcing functions 
are identical in all cases. Both kinetic and magnetic energies 
are similar in amplitudes. The biggest disparity will be noted 
in Fig. 1 10b . where the growth rates of the global mean-square 
vector potential differ significantly, resulting in A for the DNS 
run remaining about a factor of 2 larger than for any of the al- 
pha approximations at the end of the runs. 

Most of the discrepancy is accounted for by the values of 
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FIG. 13: Contour plots of the (smoothed) vector potential at t = 400 
for inverse cascades (top: DNS; middle and bottom: runs 11 and 
10 respectively). By that time, similar large-scale structures have 
formed in all runs, but note that they have different locations in the 
DNS and alpha runs. 



the k m in = 1 modes, the fundamentals, which are known 
from previous MHD inverse cascade computations to run 
away, at long times, from their nearest neighbors until lim- 
ited by their own dissipation rates l37il . Throughout the rest 
of the spectra below the forcing band, the disagreement is not 
so severe, as seen in Figs. II lfa (vector potential spectra). fTTb 
(magnetic energy spectra) and II lb (kinetic energy spectra). 
The resolutions, values of a, and plotting conventions are the 
same as in Figs. The magnetic potential is seen to follow a 
power law in agreement with the Kolmogorov-like estimation 
derived in Ref. 1 34], viz. ~ fc~ 7 / 3 ; the ensuing spectrum for 
the magnetic energy is almost flat, ~ A; -1 / 3 . The kinetic en- 
ergy spectrum on the other hand follows approximately a k 1 / 3 
law, also found in 1 26]; because the velocity field is (partially) 
slaved to the magnetic field (kinetic energy is ~ 20% of its 
magnetic counterpart), at longer times, Alfven waves can put 
in rough equipartition the kinetic and magnetic modes, except 
in the forcing band, and it is expected that Ek ~ E m below 
the forcing band. We also note that, at the onset of the inverse 
cascade process at early times, the spectra develop at large 
scales first a k 5 spectrum, followed in time by a fc 3 spectrum 
(not shown), both for the full MHD run and for the alpha runs. 
These spectra are expected because of back-scatter (see also 
Ref. 1 26] for similar results). 

Comparably wide disparities are visible in the spectral flux 
function n^/c) for the squared vector potential, shown as 
functions of wavenumber k in Fig. ^]at time t = 400. Defin- 
ing first the transfer of magnetic potential TA(k), in Fourier 
space, we have as usual for the DNS run: 

T A (k) = J At ■ JTQiJr, A})^ + c.c. , (29) 

with a similar definition for the alpha model; c.c. means com- 
plex conjugate, all fields in the above equations are taken to 
be smooth, the * indicates complex conjugate, and stands 
for taking a Fourier transform of the whole Jacobian bracket. 
The flux is then defined as usual from the transfer as: 

U A {k)= I d P T A (p). (30) 
Jo 

The plotting conventions are the same as in Figs. 1 1 01 and ITTI 
None of the alpha-model attempts comes close to the true 
MHD flux function (solid line), a phenomenon exemplified 
by the linear scale used here. 

In Figs. we display contour plots of curves of constant 
A Sz at several different times. The small-scale waviness of 
the large scale contours is due to the forcing near k ~ 20. We 
can see that the details of the A s _, -contours cease to be well 
reproduced by the alpha-model approximations, despite the 
agreements through most of the long-wavelength parts of the 
spectra. As in the earlier cases, the accurate reproduction by 
the alpha model seems to be confined to the long-wavelength 
spectral components, but does not include the detailed loca- 
tions of individual features in configuration space. 

Finally, note that the inverse cascade phenomenon in two- 
dimensional Navier-Stokes turbulence has been studied using 
the alpha model in Ref. 1 38] where it is argued that, in that 
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FIG. 14: Residual energy spectrum Ek(k) — E m (k) normalized by 
the total energy Ek(k) + E m (k) at t = 400, with plotting conven- 
tions as in Fig. 1111 The two vertical arrows correspond to the two 
values of alpha, one of which is in the middle of the forcing band 
(run 9, dashed line). Note the strong dominance of kinetic energy in 
the small scales for all runs except the DNS run (solid line) and for 
the alpha run in which no smoothing occurs for the magnetic field, 
i.e. with a m — (dash-triple dotted line). 
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TABLE II: Errors (see eq. <33» for selective decay alpha runs 14—16, 
the subscripts m and k indicating respectively, the error computed on 
the magnetic and kinetic energy spectra. 

case, the inverse cascade of energy to large scales is enhanced, 
quite substantially, when the alpha model is switched on; how- 
ever, the fact that in Ref. 1 38] both a sink at long wavelengths 
and hyperviscosity at short wavelengths are utilized may be at 
the source of this difference in behavior, compared to the case 
studied here. In that spirit, we can examine the kinetic energy 
defect in the small scales defined as 

R(k) = E k {k) - E m (k) , (31) 

normalized by the total energy density Ek(k) + E m (k) at 
that wavenumber, with Ek and E m the kinetic and magnetic 
energy spectra; as usual, when a ^ 0, the spectral density 
corresponding to the invariant (here, the total energy) will 
be taken. This energy imbalance ( I31> . when integrated over 
small scales, is known to be a source of the large-scale mag- 
netic enhancement exemplified by the inverse cascade of mag- 
netic potential, through a mechanism akin to a negative eddy 
diffusivity (here, resistivity), as shown in Ref. IB4I1 using a 
second-order closure of turbulence, or in |35]-|36] using an 



FIG. 15: Averaged (a) enstrophy, (b) square current and (c) total 
dissipation for all runs for large-scale decay (runs 13-16, see Table 
Q}. Note that the maxima are reached at the same time, but with an 
over estimation of gradients in the alpha runs. 



assumption of scale separation. In other words, it represents 
non-linear non-local interactions in Fourier space, but it does 
not say anything about the local interactions themselves. In 
Figure [T3I we plot R(k)/[Ek(k) + E m (k)] for several runs, 
with as usual the solid line for the MHD run and the dash- 
triple-dot line for the alpha run without smoothing of the mag- 
netic field. Note first that in all the runs, the fundamental mode 
is completely dominated by the magnetic energy, and in all the 
alpha-runs except the one without smoothing of the magnetic 
field (a m = 0), the small scales are completely dominated 
by the velocity; this effect is enhanced by the normalization, 
but little energy resides in the smallest scales, beyond alpha, 
so that the inverse cascade can still take place in alpha runs, 
albeit at a slower pace. In other words, we see that in the 
framework of the alpha model, this energy defect is modi- 
fied from the usual MHD case, and more importantly it even 
changes sign for the larger values of a at high wavenumber, 
becoming positive and thus indicative of too high a dissipation 
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FIG. 16: Averaged (a) magnetic energy and (b) kinetic energy spectra 
for the large-scale turbulence runs 13-16 (see Tabled between t = 3 
and t = 7 (first and second peak of the enstrophy). The wavenumber 
corresponding to a is indicated by the vertical line, and the under- 
resolved run is shown with a dotted line. At these early times, spectra 
are in close agreement except at small scales. 

of the magnetic field in the presence of alpha-smoothing of 
small scales, i.e. with a m — 0. Even though alpha modeling 
is about the dynamical evolution of scales larger than alpha, 
the scales smaller than alpha nevertheless participate into the 
dynamics of the flow including at large scale and are responsi- 
ble for the discrepancy in the growth rate of squared magnetic 



potential we observe here (see FigllOb). Thus, such an effect 
presumably linked to the fact that the magnetic field is decay- 
ing with a scale dependency that is leaning more heavily on 
the small scales (like an effective fc 4 -type hyper-resistivity), 
may be at the source of the different behavior we observe be- 
tween the different alpha-runs and the DNS run; if only non- 
local transfer (in Fourier space) is active in the inverse cas- 
cade in the case of the alpha runs, this could be at the origin 
of the slower growth of A. This point will require further in- 
vestigation in the three-dimensional case, in relation with the 
large-scale dynamo problem, since the growth rate associated 
with the inverse cascade of magnetic helicity is the relative 
(kinetic helicity minus magnetic current helicity) in the small 
scales [29]. In the limit of very large a, the dynamics become 
trivially simple; see Ilyin and Titi |39|. 



Two separate questions of accuracy are addressed in this 
Section. First, we attempt quantitative comparisons of alpha 
model computations run on a 256 2 grid with the standard pro- 
vided by a well-resolved MHD run at a resolution of 1024 2 
for the same initial conditions and the same Reynolds num- 
bers. Secondly, we compare the same alpha-model runs with 
DNS solutions on a 256 2 grid in which no smoothing has been 
applied and which are undoubtedly unresolved, in that the dis- 
sipation wavenumbers (estimated on the basis of viscous and 
Ohmic enstrophy dissipation) exceed the maximum values of 
k retained in the computation; such runs correspond to under- 
resolved computations. This latter test is necessary if one is 
to argue that accuracy has been improved by the alpha model 
in computations of comparable resolution. Error estimates are 
provided first for freely-decaying turbulence at early times, 
near the peaks of the dissipations (runs 13-16 of Table[Q, and 
then for the situations described in Sees. IV and V, where later 
times (and hence greater accumulated errors) are involved. 

First we examine the case of freely decaying turbulence 
at an early time, which could be considered to be the early 
stages of a selective decay computation. The turbulence ini- 
tially is confined to a band of wavenumbers between k% = 1 
and k,2 — 3. In Figure ^] we show the computed mean- 
square vorticity (a), mean square current density (b), and mean 
square dissipation v < oj 2 > +r/ < j 2 > as functions of 
time. We show runs under all four circumstances, starting 
from the same initial conditions: fully-resolved 1024 2 MHD, 
an alpha model run at 256 2 with a = 1 /50 for both fields, an 
alpha model run with an unsmoothed magnetic field (a m = 0) 
and «fc = 1/50 for the velocity field at 256 2 , and an under- 
resolved but unsmoothed DNS run at 256 2 ; we see that time 
scales of growth of field gradients, and the amplitudes of such 
gradients are comparable but not identical. 
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TABLE III: Errors as defined in eq. l33lfor the inverse cascade alpha runs 9-12 (see Table[I}; E m and E k indicate respectively the magnetic 
and energy spectra. Note the larger errors for run 9 corresponding to the alpha wavenumber (k a — 20) embedded in the forcing band. 



Note that for the fully resolved DNS run, k rnax = 341 to 
be contrasted with the dissipation wavenumber based on the 
magnetic variables (computed using a Kolmogorov spectrum) 
of kdiss ~ 280 at the peak of dissipation i.e. for t ~ 7 (and 
kdiss ~ 250 when based on the velocity); this shows that 
the run is well resolved at all times, and indeed no "bottle- 
neck" appears in the compensated spectra, where by bottle- 
neck is meant an accumulation of energy in the small scales, 
seen as a bump in the compensated spectra for k ~ k max ; 
this phenomenon has sometimes appeared in the presence of 
under-resolved computations, or when using hyper-viscous or 
hyper-resistive dissipative operators 12611 . At the same time 
t ~ 7, the Taylor Reynolds number (see equation (I25H is 
R\ ~ 1150, with a Taylor wavenumber k\ ~ 20 (with 
k\ ~ 17 when based solely on the current density); in this 
case, only run 14 (with both alphas non zero) has a slightly 
higher R\ at that time, whereas runs 15 (with a m = 0) and 
run 16 (under-resolved DNS) have slightly lower R\. 

Figure[H)]shows the energy spectra, for the freely decaying 
run, averaged over the period of time when the dissipation 
remains quasi-constant (t = 3 to t = 7), and Fig. I17ldisplavs 
the relative error spectrum, defined as 



e r (k) 



E*(k)-E(k) 

W) 



(32) 



Here, E(k) is the well-resolved energy spectrum (either ki- 
netic or magnetic), summed over all /c-vectors with the same 
magnitude, regarded as the "truth," and E* (k) is the energy 
computed from any one of the stated three approximations to 
it. 

Both figures show that the under-resolved flow displays sys- 
tematically a greater error than the alpha model run, at large 
scales (k = 1 to k = 3) as well as at small scales near the 
cut-off, 1/a. Finally, in Fig. [^] are plotted the pdfs of (a) the 
current and (b) the spatial density of the dissipation of energy, 
the former in lin-log scales, and the latter in log-log scales. 
The same line conventions as in Fig. [H)]have been adopted. 
The alpha model clearly reproduces better the gradients, at a 
given resolution, than the under-resolved flows. 

Two possible global quantifications for measuring the er- 
rors are El and E2, defined by the relations: 
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E(k) 
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Thus a small value of E\ or E2 will indicate a closeness on 
the part of the MHD approximations (alpha-modeled or unre- 
solved) to the full MHD DNS results. 

Ei and E2 are exhibited as TablellTlfor the freely decaying 
runs just described (and indicated by the superscript "D"). The 
conventions used in the Table are that the number of the run is 
followed by a subscript, where m stands for magnetic spectra 
whose error is being assessed, and k for kinetic ones. Time 
averages have been performed from t = 3 to t = 7, the vicin- 
ity of the main peaks in the dissipation rates. This is thought 
to be the time when the turbulence is of its most broad-band 
character, when the alpha model would be having its strongest 
impact. 

The errors so defined are smaller for these early times for 
this freely decaying situation than they are for the situations 
studied in Sections III, IV, and V which examine late-time 
evolutions. We note that the lowest error (by a factor of ~ 2 
compared to the unresolved run) occurs for the alpha run with 
both alphas equal (run 14, see Table [j}. Note also that Fig. 
[n] shows errors at each k, whereas the errors in Table HP are 
normalized by the total energy (truncated at 1/a) and hence 
are smaller. 

In Table IIHI we show the errors E\ and E2 for the inverse 
cascade situation of Sec. V and as indicated by the super- 
script I; the runs are the same as those in Sec. V, and given 
by their number (see Table []}, with E m and Ek standing as 
usual for magnetic and kinetic energy spectra. These errors 
are computed at late times, when the true solution may be ex- 
pected to have drifted further from the 1024 2 run and hence 
lead to larger errors than in the freely decaying runs of Table 
UTI Moreover, as expected, when the alpha cut-off is too close 
to the forcing band, the errors are larger, both for the kinetic 
and the magnetic spectra. 

Finally, in Table HV1 the E\ and E2 errors are shown for 
the dynamic alignment runs 6 and 7 of Sec. IV, with a super- 
script "DA" and a supplementary index "e" or "1", in order to 
indicate early or late times in the run; specifically, early signi- 
fies that the average is taken for for 5 < t < 10, and late for 
55 < t < 60. Again, the type of spectrum for which the error 
is displayed (either E m or Ek) is given before the run number. 
The low-resolution computation (run 7) does not have signif- 
icantly higher errors at early times, but errors accumulate at 
later times, more so for the lower resolution computation (run 
7). The reason for time averaging the magnitude of the er- 
rors is that, when plotted as functions of time, the error curves 
cross each other repeatedly. Time intervals can be found when 
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either one is smaller than another. Time averaging, over inter- 
vals long enough to contain many of these crossings but short 
compared to the duration of the runs, has seemed to provide 
the most objective number for addressing which error is "typ- 
ically" smaller. 



VII. VERY HIGH REYNOLDS NUMBERS 

The eventual utility of the alpha model if it can be justified 
will be that it will permit explorations of Reynolds number 
regimes that are far above those that can be obtained from di- 
rect numerical solutions of the MHD equations. It may be 
noted that extimates have been given for the number of de- 
grees of freedom of the Navier-Stokes alpha model lEoll : see 
also fT2ll . Whereas most of this paper has been devoted to 
regimes in which direct MHD solutions can be compared to 
alpha model solutions, we have thought it interesting to show 
one alpha model computation that goes beyond what can be 
contemplated from unsmoothed solutions presently. We do 
this without any definitive claims for accuracy, but just as a 
suggestion of what the alpha model might provide in the way 
of future predictions. A detailed analysis of high Reynolds 
number runs at higher resolutions than what is performed here 
and using the alpha model will be presented elsewhere. 

We display results for a run with 77 = v = 2 x 10~ 5 , a time 
step of At = 2.5 x 1CT 4 , and 2048 2 grid points. The value 
of a -1 is chosen to be 300 ~ k max /2, with k max ~ 682 
for this run. The initial (equal) kinetic and magnetic energies 
are loaded with random phases into the ring from k\ = 1 
to k 2 = 3 in Fourier space, and the r.m.s. values of v and 
B are unity. Computed at t = 7, i.e. close to the maximum 
of dissipation (see Fig. 119b). the Taylor Reynolds number 
R\ ~ 5200. This is roughly comparable with what can be 
accomplished with a DNS on a grid of more than 10 8 points. 

Fig. 119b shows the evolution of the total kinetic energy 
(dashed line) and total magnetic energy (solid line), referred 
as usual to unit volume, as functions of time. Fig. 119b shows 
the evolution of the mean square vorticity (dashed line) and 
mean square current density (solid line) as functions of time. 
The qualitative behavior of all quantities in Figs. ^]will be 
seen as not significantly different from that observed at lower 
resolutions, although oscillations in the kinetic and magentic 
energies are persistent until t ~ 6. But the idea has been to 
extend the inertial range as much as possible. It will be noted 
that the magnetic and kinetic energies are far from equiparti- 
tioned, with magnetic energy in excess by a factor ~ 3 by the 
end of the run. The peak value of mean dissipation at t ~ 7 
(e ~ 0.074) is comparable with that at the lower Reynolds 
number (runs 13-16), confirming previous results (see e.g. Fig. 
7 in Ref. l4lll ) of lack of dependency of e with Reynolds num- 
ber, at least at a magnetic Prandtl number of unity. 

Compensated energy spectra are displayed in Figs. |20] Fig. 
I20h is the kinetic energy spectrum, multiplied by fc 5 / 3 , and 
averaged between t = 2 and t = 6, i.e. the times at which 
the energies are oscillating with little dissipation yet. Fig 
120b shows the magnetic energy spectrum averaged over the 
same time interval and similarly multiplied by fc 5 / 3 . Fig. 120b 
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FIG. 17: Normalized errors, shown up to a -1 , in time-averaged (a) 
magnetic and (b) kinetic energy spectra for the large-scale turbulence 
runs 13-16. Larger discrepancies occur both at small and at large k, 
in particular for the unresolved run (dotted line). 
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TABLE IV: Errors on the magnetic E m and kinetic Et energy spec- 
tra for dynamic alignment runs (see Table|3 for early (e) and late (1) 
phases of the evolution. As expected, errors are larger at later times, 
and are larger for the low resolution computation (run 7). 



shows the total energy spectrum, similarly compensated and 
time averaged over the same time interval. In all three fig- 
ures, the horizontal dotted line has zero slope, and so would 
coincide with a fc~ 5 / 3 Kolmogorov inertial range spectrum so 
compensated. In all three figures, the dashed line has slope 
1/6, and so would be tangent to a A; -3 / 2 spectrum (as pro- 
posed by Iroshnikov and Kraichnan 1 42, 43]) which had been 
multiplied by fc 5 / 3 . It would appear from these figures that 
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/c~ 5 / 3 would fit the computed spectra significantly better than 
/c~ 3 / 2 would. We do not attach any finality or conclusiveness 
to this observation, because intermittency is known to steepen 
energy spectra obtained from dimensional analysis, and it is 
known that high-order structure functions computed for sta- 
tistically steady 2D-MHD flows display a behavior that dif- 
fers from that of turbulent neutral (3D) fluids 1 44-] . The total 
energy spectrum computed with the alpha model agrees with 
other findings (see e.g. 1 44, 45] [26]) at lower Taylor Reynolds 
numbers. In that light, we conclude that the alpha model does 
not alter previously known results, suggesting that intermit- 
tency is worth further investigation in the context of the alpha 
model, both in two and three dimensions. In the same spirit, 
we note the absence of any "bottle-neck" in the spectra, even 
at these high Reynolds numbers, in contrast to what is found 
in jH. 

Finally, it is worth noting that, computed at the maximum 
of dissipation at t ~ 7, the Taylor wavenumber in the al- 
pha run is k\ ~ 55, i.e. well below the alpha wavenum- 
ber of 300, whereas the dissipation wavenumber, based on 
the alpha-enstrophy < lolo s > and the square current < 
j 2 >, are respectively 1300 and 1500, i.e. well beyond the 
largest resolved wavenumber (k max = 682) for this compu- 
tation; however, even when plotting the dissipation spectrum 
k 2 [Ek(k) + E m (k)], no bottle-neck is observed (not shown). 



VIII. SUMMARY AND DISCUSSION 

The intent of this article has been an empirical study of 
the extent to which the alpha model, or Lagrangian-averaged 
model, equations predict the results of computed incompress- 
ible 2D MHD behavior in which no further modeling or ap- 
proximations are made. Our primary motivation has been to 
acquire confidence in the alpha model in hopes that it can be 
used for problems with such high Reynolds-like numbers that 
they cannot be computed in the framework of the primitive 
MHD equations. A useful starting point has seemed to be 
to work on classical problems in rectangular periodic condi- 
tions where some information about solutions has accumu- 
lated over the last thirty years. These problems include those 
often grouped under the terms selective decay, dynamic align- 
ment, direct and inverse cascades, and tabulating frequency 
distributions or pdfs for fluctuating field quantities. We have 
compared the results of direct numerical solutions for these 
problems with the solutions of alpha-model equations for the 
same initial and boundary conditions, both for freely decay- 
ing and for forced turbulence. In order to do so at sufficiently 
high Reynolds numbers, we have kept our comparisons to the 
two-dimensional geometry for which reasonable resolutions 
can be obtained without overly taxing available computer re- 
sources. 

The principal success we have to report is in the satisfac- 
tory reproduction of the long-wavelength spectral behavior 
for magnetic and velocity field evolution, in nearly all cases. 
Also, characteristic length scales of the flow, as the Taylor 
scale, are well reproduced by the alpha modeled simulations. 
The only exception is that some discrepancies remain for the 
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FIG. 18: Normalized pdfs of (a) the current density, and (b) ln(^o; 2 + 
r)j 2 ), for the large-scale turbulence runs (see Table[]}. At this time, 
small-scale structures are less developed in the alpha runs than in the 
DNS run (solid line). 



longest-wavelength behavior in the case of driven cascades 
of mean-square magnetic vector potential at late times. The 
model has not proved accurate in the detailed reproduction of 
specific spatial features of the turbulence at earlier times. Nor 
does the alpha model reproduce accurately the pdfs of inter- 
mittent fluctuations of large amplitude, and likely because of 
the deliberate suppression of small spatial scales, is not ex- 
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FIG. 19: (a) Magnetic and kinetic energy, and (b) enstrophy and 
square current as a function of time for run 17 (see Table []} with 
2048 2 grid points and a -1 = 300. Note the several peaks in the 
enstrophy and the mean square current density, as well as the excess 
in magnetic excitation, both at large scale and at small scale. 






FIG. 20: (a) Kinetic energy spectrum, (b) magnetic energy spectrum, 
and (c) total energy spectrum compensated with fc 5 ^ 3 , for run 17 (see 
Table[]}. Temporal average is performed between t = 2 and t = 6. 
The horizontal dotted line corresponds to a fc _5//3 spectrum, and the 
dashed line to a k~ 3 ' 2 law. A clear inertial range extends for more 
than one decade in wavenumber. 



pected to. 

It should be remarked that there is the more difficult prob- 
lem of obtaining a clear physical understanding of the nature 
of the alpha approximation itself. It was originally arrived at 
by mathematics of considerable sophistication and complex- 
ity, which left intuitive gaps in just what was being assumed. 
This is a very different perspective than the one used in Ref. 
lll3ll or Section II of this paper, where the recipe of smoothing 
the v and B fields but not their sources, and then neglect- 
ing the fluctuations in those fields about their averages was 
invoked. Neither prescription seems clear enough to us at a 
physical level to argue for it strenuously on any basis other 
than its satisfactory consequences. Nor is it clear why two 
such apparently different procedures should end at the same 
place, despite the happy fact that they seem to do so. Clari- 



fying the conceptual foundations of such modeling at an intu- 
itive physical level is a serious but stimulating challenge. 

In summary, our judgment of the alpha model is that it re- 
produces satisfactorily the time development given by well- 
resolved DNS computations for spectra up to about k ~ 1/a 
and does well enough at reproducing the probability distribu- 
tion functions of fluctuations. What it does not do satisfacto- 
rily is to reproduce the locations, trajectories, and shapes of 
structures in configuration space. 

We see several directions in which to extend this work. For 
example, there is the natural one of three-dimensional com- 
putations, still in rectangular periodic boundary conditions, 
where such problems as the inverse cascade of magnetic helic- 
ity (an inherent part of the large-scale dynamo problem), the 
spectral anisotropy induced by the presence of a dc magnetic 
field 1 46, 47], and the small magnetic Prandtl regime remain 
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to be investigated (see for recent studies in the latter case e.g. 
1481 149[). Finally, the questions of material boundaries with 
non-ideal boundary conditions and departure from rectangular 
to spherical or cylindrical symmetry seem necessary as well. 
The only work published so far involving material boundaries 
and the alpha model seems to be that of Chen et al. 0.0]. It 
is our intent to move in these directions in the near future. 
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